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Abstract 

A non-singlet QCD analysis of the structure function F2(x,q 2 ) in LO, NLO, NNLO and 
N 3 LO is performed based on the world data for charged lepton scattering. We determine 
the valence quark parton densities and present their parameterization and that of the cor- 
related errors in a wide range of x and Q 2 . In the analysis we determined the QCD-scale 
A QCD,N f =4 = 265±27MeV (NLO), 226±25MeV (NNLO), 234±26MeV (N 3 LO),witha 
remainder uncertainty of ±2 MeV for the yet unknown 4-loop anomalous dimension, corre- 
sponding to a s (Mf) = 0.1148±0.0019 NLO, 0.1134 IS NNLO, 0.1141 ^022 N 3 LO. 
A comparison is performed to other determinations of the QCD scale and a s (M^) in deeply 
inelastic scattering. The higher twist contributions of F^XjQ 2 ) and F^x^Q 2 ) are ex- 
tracted in the large x region, subtracting the twist-2 contributions obtained in the NLO, 
NNLO and N 3 LO analysis. 
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1 Introduction 



Deeply inelastic scattering (DIS) at large 4-momentum transfer Q 2 = —q 2 offers the possibility 
to probe the partonic sub-structure of the nucleon. The light-cone expansion [f] allows to 
describe the forward Compton amplitude in terms of the twist-expansion. For the twist-2 
contributions the factorization theorems hold and the calculation of the scaling violations of 
the parton distribution functions and structure functions is possible within the framework of 
perturbative Quantum Chromodynamics (QCD). In the present paper we determine the flavor 
non-singlet parton distribution functions xu v (x,Q 2 ) and xd v (x,Q 2 ) using the available e(fi)p 
and e(n)d world data [2-6] up to the next-to-next-to leading (NNLO) and next-to-next-to-next 
to leading order level (N 3 LO). This is possible due to the advent of the 3-loop non-singlet 
anomalous dimensions [7] and 3-loop Wilson coefficients [8] 2 , as well as first informations on the 
4-loop non-singlet anomalous dimension, with its second moment [10]. The QCD analysis of 
the world non-singlet data has the advantage to deal with a restricted set of parton distribution 
functions, to which the virtual photon couples directly and which do therefore contribute already 
at leading order. The analysis is, furthermore, free of assumptions on the gluon distribution 
function which stabilizes the systematics since a lower number of parameters describing the 
parton distributions has to be measured and the analysis is free of the correlation between the 
QCD scale Aqcd and the gluon density. Also due to the particular behaviour of the flavor 
non-singlet parton densities in the smaller Bjorken-x region in comparison to the sea-quark 
and gluon densities a separate non-singlet analysis is in order to better understand systematics 
effects present in combined non-singlet and singlet analyses. 

Previous 3-loop QCD analyses were mainly performed as combined singlet and non-singlet 
analyses [11, 12] 3 , partly based on preliminary, approximative expression of the 3-loop splitting 
functions. Other analyses were carried out for fixed moments only in the singlet [14, 15] and 
non-singlet case analyzing neutrino data [14] 4 . First results of our non-singlet analysis were 
published in [17]. Very recently a 3-loop non-singlet analysis was also carried out in Ref. [18]. 
There are numerous NLO analyses, among them [2,5,6,19], which we will use in comparisons 
below. Another recent NLO parameterization [20] fixes a s (M§) and presents a series of ten 
parton distribution sets associated to equidistant values of the strong coupling constant in the 
range a s {M 2 z ) e [0.110, 0.128]. 

The non-singlet analysis will provide us with the distribution functions xu v (x, Q 2 ) and 
xd v (x, Q 2 ) with correlated errors. From these distributions we may form Mellin-moments which 
can be compared to lattice-simulations of these quantities. The lattice calculations have to be 
performed at as low as possible pion masses due to a potential non-linear behaviour on m n in 
the small mass region and non-perturbative renormalization has to be used. Precision results 
of these simulations with Nf = 2 dynamical quarks are expected soon. The comparison of the 
Mellin moments obtained in a 3(4)-loop QCD analysis of the non-singlet World Data at the 
one side and lattice simulations on the other side will lead to new important tests of Quantum 
Chromodynamics. 

The paper is organized as follows. In Section 2 we briefly summarize the basic formalism. 
The data analysis and calculation of errors is outlined in Section 3. Different corrections such as 
the account for a finite sea quark difference x(d — u)(x,Q 2 ), heavy quark contributions, target 
mass and deuteron wave function corrections are dealt with in Section 4. The fit results for 

2 Moments of the 3-loop anomalous dimension and Wilson coefficients were calculated in [9]. 
3 A recent update of the latter analysis has been given in [13] very recently. 
4 For earlier results see [16]. 
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the non-singlet parton distribution functions, their correlated errors and evolution are presented 
in Section 5, where also comparisons to other analyses are made. In Section 6 the results on 
Aqcd and a s (M|) in NLO, NNLO and N 3 LO are given and compared to other NLO and NNLO 
results. In Section 7 we present our results on the higher twist contributions in the large x region 
for the structure functions F$ (x, Q 2 ) and F^(x, Q 2 ) and Section 8 contains the conclusions. 

2 Basic Formalism 

We choose the following parameterization for the valence quark densities 

xvi{x,Q 2 ) = A iX a *{l -x) bl {l +piVx + j iX ) , (1) 



with Vi = u v ,d v , i = u,d. Three valence-projections are considered. In the region x < 0.3 for 
the difference of the proton and deuteron data we use 



x(u + u) (x, Q 2 ) — x(d+ d) (x, Q 2 



(2) 



The combinations of parton densities in the valence region x > 0.3 for ' d (x, Q 2 ) are 



f;(x, Q 2 ) = — [x(u - u)(x, Q 2 ) + x(d - d)(x, Q 2 )\ + - [x(u - u)(x, Q 2 ) - x(d - d)(x, Q 2 ) 



f v A*,Q 2 ) 
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x(u — u) (x, Q 2 ) + x(d — d) (x, Q 2 ) 



(3) 
(4) 



The evolution equations are solved in Mellin-iV space and the Mellin transforms of the above 
distributions are denoted by f NS (N,Q 2 ), fp d (N,Q 2 ), respectively. The non-singlet structure 
functions are given by 



F k (N,Q 2 ) = l + a s (g^)C 1 (7V) + af(g 2 )C 2 (7V) + af(Q^)C 3 (7V) f k (N, Q 1 ) 



(5) 



for the three cases above. Here a s (Q 2 ) = a s (Q 2 ) /(47r) denotes the strong coupling constant and 
Ci(N)(Q 2 ) are the non-singlet Wilson coefficients in 0(a l s ) [8,21,22]. For the representation of 
the Wilson coefficients we used those given in [8,23] which are based on numerical fits. For the 2- 
loop Wilson coefficients and 3-loop anomalous dimensions compact analytic representations are 
also available, based on the reduction of the respective harmonic sums to 5 (14) basic functions 
[24-28]. 5 

The solution of the non-singlet evolution equation for the parton densities to 4-loop order 
reads 



F k (N,Q 2 ) = F k (N,Q 2 Q )(- 

\a 



a \ -A>(A0//3o 



l-^(a-ao) 



Po 



5 Similar representations were derived for the 2-loop Wilson coefficients for the unpolarized and polarized 
Drcll-Yan process and scalar and pseudoscalar Higgs-boson production at hadron colliders in the heavy mass 
limit as well as for the hadron fragmentation functions in e + e~ annihilation [29,30]. 
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Here, Pk denote the Mellin transforms of the (k + l)-loop splitting functions. 
The strong coupling constant obeys the evolution equation 



da s (Q 2 
dlnQ 5 



Y.Pk a 

k=0 



k+2 

s 



(Q 2 ) ■ 



(7) 



Eq. (7) is solved in the MS-scheme applying the matching of flavor thresholds at Q 2 = m 2 c and 
Q 2 = ml with m c = 1.5 GeV and = 4.5 GeV as described in [31,32]. The convention for the 
MS-scheme introduced in [33] is extended in this way. To be capable to compare with other 
measurements of Aq CD we adopt this prescription. 



3 Data Analysis and Error Calculation 

For the non-singlet QCD analysis presented in the present paper structure function data mea- 
sured in charged lepton proton and deuteron deep inelastic scattering are used. The experiments 
contributing to the statistics are BCDMS [2], SLAC [3], NMC [4], HI [5], and ZEUS [6]. The 
BCDMS data were recalculated replacing Rqcd by -R1998 [34]. The SLAC data come from a 
re-analysis of the data from the experiments E49, E61, E87, and E89. All deuteron data were 
corrected for Fermi motion and off-shell effects [35]. We form three data samples : F^x^ 2 ), 
F^ix, Q 2 ) in the valence quark region x > 0.3 and F 2 NS = 2(F 2 — F 2 ) in the region x < 0.3. 
In the valence quark region we approximate the parton distribution functions by pure valence 
quarks. Only data with Q 2 > 4 GeV 2 were included in the analysis and a cut in the hadronic 
mass of W 2 > 12.5 GeV 2 was applied in order to widely eliminate higher twist effects from the 
data samples. After these cuts we are left with 762 data points, 322 for F 2) 232 for F 2 , and 
208 for F 2 NS . The BCDMS data have been used with an additional cut of > 0.3 in order to 
exclude a region with significant correlated systematic errors, cf. also [5]. In case of e(n)N deep 
inelastic scattering experiments measuring the kinematic variables either only form the lepton 
or the outgoing jet-kinematics the systematics of the respective measurement becomes worse 
at low values of y (lepton measurement) or high values of y (jet measurement). These effects 
were studied in Ref. [36] and give rise to the y-cut mentioned for the BCDMS-data. A cut of 
Q 2 > 8 GeV 2 was imposed on the NMC data. These cuts reduce the total number of data points 
available for the analysis from 762 to 551. 
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The data sets used contain both statistical and systematic errors which in turn, for some data 
sets, split into uncorrelated and correlated ones. There are several procedures proposed to treat 
the experimental errors and the determination of uncertainties from parton distributions [37]. 
Here we decided to use the simplest procedure by adding the statistical and systematic errors, 
uncorrelated and correlated ones, in quadrature when using the data sets in the fit. However, we 
allowed for a relative normalization shift iVj between the different data sets within the normaliza- 
tion uncertainties AiVj quoted by the experiments or assumed accordingly. These normalization 
shifts were fitted once and then kept fixed. The number of data points, their x and Q 2 range, 
and the normalization shifts determined are summarized in Table 1. 

The normalization shift for each data set enters as an additional term in the x 2 -expression 
which then reads 

" ( N t - l) 2 (NiFff - Fi h f or ' 2 



x 2 = E 



+ E 



(8) 



(AiV,) 2 fri (N t AF^ ta y 

where the sums run over all data sets and in each data set over all data points. The minimization 
of the x 2 value above to determine the best parameterization of the unpolarized parton distri- 
butions is done using the program MINUIT [38]. Only fits giving a positive definite covariance 
matrix at the end have been accepted in order to be able to calculate the fully correlated la 
error bands. 

The la error for the parton density f q as given by Gaussian error propagation is 

where the sum runs over all fitted parameters. The functions df q /dpi are the derivatives of 
f q w.r.t. the fit parameter pi, and cov(pi,pj) are the elements of the covariance matrix. The 
derivatives df q /dpi can be calculated analytically at the input scale Q\. Their values at Q 2 are 
given by evolution which in our case is performed in Mellin-iV space. The evolution formalism in 
Mellin-iV space is such that the input and the evolution parts factorize which leads to a straight- 
forward analytic error calculation. The derivatives evolved in Mellin-iV space are transformed 
back to x-space and can then be used according to the formula above together with the elements 
of a positive definite covariance matrix determined at the input scale Qq by the QCD fit to the 
data. 

Let us discuss the derivatives in Mellin-iV space a bit further. The Mellin-iV moment for 
complex values of N calculated at the input scale Ql for the parton density parameterized as in 
Eq. (1) is given by 



M\f(x,a,b,p,i)\(N) = A[B(a + N + 1,6 +l)+ P B(a + N -1/2,6+1) 

+ -fB(a + N, b + 1)] , (10) 

with the normalization constant 

C 

A — (id 

[B(a, b+l)+ P B(a + 1/2, b + 1) + 7 5(a + 1, b + 1)] ' y ! 

Here C is the respective number of valence quarks, i.e. C Uv = 2 and Cd v = 1, and B(a,b) = 
r(a)r(6)/r(a + 6) denotes the Euler Beta-function for complex arguments. 
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The general form of the derivative of the Mellin moment M w.r.t. the parameter p is given 



by 



dM[f(x,p)](N) dM OA 

Op Op Op 



(12) 



with M = M/A. In the present analysis only the parameters a and 6 have been fitted for both 
the u v and the d v parton distribution while the other parameters involved are kept fixed after a 
first minimization since their errors turned out to be rather large compared to the central values. 
For the parton densities u v (x, Ql) and d v (x,Ql) the derivatives w.r.t. parameter a are: 



dM 

da 



{[^( a _ i + N ) _ + N + 6 )] + 7 ^_i±_^ (^( Q + 7V)- 



*(a + N + b + l))}B(a - 1 + AT, b + 1) + p[W (a - - + iV 
* ( a + 2 + N + b )^ B (° ~ \ + N ' b + 1 



and 



dA 

— = -AZ a /X a = -CZ a /X 2 a , 



where 



and 



B(a, b + 1) + pB (a + ^ b + lj + 7 5(a + 1,6+1). 



For the derivatives of u v (x, Q%) and d v (x, Q 2 ,) w.r.t. parameter b one obtains: 



dM 
~db 



{[^(6+l)-^(a + iV + 6)]+ 7 



a-l + N 



(*(&+!)- 



a + N + b 

V{a + N + b+ l))}B{a - 1 + N, b + 1) + p [*(6 + 1)- 
* (a + i + iV + bj 5 (a - ^ + iV, 6 + l) 



and 



(13) 



(14) 



[*(a) - (a + 6 + l)]5(a, 6 + 1) + p * (a + ^ - * (a + ^ + 6 + l) 
S(a + -, 6 + 1) + 7[*(a + 1) - *(a + 1 + 6 + l)]B(a + 1,6+1), (15) 



(16) 



;i7) 



(9/1 

— = -AZ b /X b = -CZ b /X 2 b , 



(18) 
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where 



Z h = 




(19) 



and 




(20) 



Here ^(z) = dlnT(z)/dz is Euler's ^-function. To obtain the gradients for the error calculation 
of the structure functions and _F 2 NS the relevant gradients of the parton distribution 

functions in Mellin space have to be multiplied with the corresponding Wilson coefficients. 

This yields the errors as far as the QCD parameter A QCD is fixed and regarded uncorrelated. 
The error calculation for a variable Aq CD is done numerically due to the non-linear relation and 
required iterative treatment in the calculation of a s (Q 2 , Aq CD ). The respective gradients are 
given by performing the evolution for A ± 5, with J < A, and calculating 



using typical values for 5 < 10 MeV in the present analysis. 

4 Corrections 

To have the possibility to include deuteron data into the analysis deuteron wave function correc- 
tions both for Fermi motion and off-shell effects [35] have to be carried out. These corrections 
are performed during the fit. 

Since we analyze charged lepton-nucleon scattering the non-singlet structure functions re- 
quire the distribution x(d — u)(x, Q 2 ) as an input, which cannot easily be extracted in the DIS 
analysis, but can be determined from Drell-Yan data [39] . We refer to a parameterization in [40] 
at Q 2 = 1 GeV 2 which has been evolved to Ql = 4GeV 2 . An illustration is given in Figure 1, 
containing also the older parton parameterization [12] which has been updated recently [13]. 

Part of the data is situated at low values of Q 2 and large values of x which have to be 
corrected for target mass effects, cf. [41]. The respective numerical effects are discussed below, 
see Figures 5, 6 and 8. 

Beginning with 0(a 2 ) the non-singlet structure functions receive heavy flavor contributions 
due to the Compton-process 7* + q(q) — > q(q) + <?*—>■ QQ and 7* + g* — > QQ, where the 
off-shell gluon is emitted from a light quark or anti-quark line in lowest order, [42]. As the 
present analysis is carried out in Mellin space one needs accurate and fast representations for 
the heavy flavor contributions, which were derived in Ref. [43] for all available heavy flavor 
Wilson coefficients. We accounted for these contributions in the present analysis. The numerical 
effects of these contributions for the case of cc production are shown in Figures 2, 3 and 4 for the 
proton and deuteron valence regions x > 0.3 and for -F 2 NS (x, Q 2 ) at x < 0.3 for a wide region in 
Q 2 . The respective arrows mark the lowest data point in the analysis. Except of the very high 
Q 2 region, Q 2 ~ 10000 GeV 2 , where the relative contribution amounts to O(0.4%) the typical 




(21) 
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corrections are (far) below 1 per mille and are therefore widely negligible. The present analysis 
would require the 0(af) heavy flavor corrections for -F 2 NS (:r, Q 2 ) as well, which were not yet 
calculated. 6 From the 0(a 2 ) results we do not expect large additional effects in the non-singlet 
case, however. The 66-contributions are even smaller. 

Universal small-rr contribution to the evolution kernel for deep-inelastic scattering were de- 
rived in [46] in leading order. Resumming these effects with renormalization group techniques 
yields contributions on the level of 1-2% and less depending on the size of the subleading 
terms [47], which are important. Like in the flavor singlet case [48] the sub-, subsub-leading 
contributions, etc. turn out to be of quite similar size, however, with alternating sign, as the 
respective next term. This is easily verified expanding the available fixed oder results in con- 
secutive powers in Mellin 1/N k , cf. also [7]. Hence, these effects do not affect the present 
analysis. 

5 Fit Results 

The flavor non-singlet analysis relies on the three complementary data sets: the structure func- 
tions F2 ,d (x, Q 2 ) in the valence-quark region x > 0.3 and the combination of these structure func- 
tions F 2 NS (:r,Q 2 ) = 2[F 2 D (x,Q 2 ) - F$(x,Q 2 )] . To isolate the twist-2 contributions to the deep- 
inelastic structure functions F2 ,d (x, Q 2 ) cut-studies in Q 2 and W 2 were performed. It turns out, 
that power corrections are widely absent in the kinematic region Q 2 > 4 GeV 2 , W 2 > 12.5 GeV 2 . 
We therefore choose these cuts for the whole data set to perform a twist-2 non-singlet QCD 
analysis. We further unfold the target mass corrections [41] for all data. 

In Figure 5 the proton data for F2(x,Q 2 ) are shown in the valence quark region x > 0.3 
indicating the above cuts by an arrow. The solid lines correspond to the NNLO QCD fit and the 
dashed lines correspond to the NNLO QCD fit adding target mass corrections. We extrapolate 
these fits to the region 12.5 GeV 2 > W 2 > 4 GeV 2 . There, at higher values of clear gap 
between the data and the QCD fit is seen. Figure 6 shows the corresponding results for the 
deuteron data. In the kinematic region x < 0.3 the above W 2 cut is ineffective. Figure 7 shows 
the result of the NNLO fit for this region. The effect of target mass corrections in the latter 
kinematic domain is shown in Figure 8, which turns out to be rather small. 

In Table 2 we summarize the fit results for the parameters of the parton densities xu v (x, Qq), 
xd v (x, Qq) and Aqcd i n the MS-scheme in NLO, NNLO, and N 3 LO. The value of x 2 / n df a t the 
minimum improves from 0.89 at NLO to 0.84 at N 3 LO. In all the fits the covariance matrix is 
positive definite. As an example we show the covariance matrix of the NNLO fit in Table 3. 

Figure 9 illustrates our fit results for xu v (x,Ql), xd v (x,Ql) at Ql = 4 GeV 2 at NNLO with 
correlated errors (left figures). We compare to the results of [49] and a very recent analysis [13]. 
While there is very good agreement in the case of xu v (x, Qq) some differences at the la level are 
visible for xd v (x, Ql)- We also compare our results to some NLO analyses [5,19,50] (right figures). 
The spread here turns out to be somewhat larger than in the NNLO case. Note that the analyses 
we compare to are combined non-singlet/ singlet analyses and already different assumptions on 
the sea quark distributions may cause deviations for the non-singlet distributions in the range 
of smaller values of x. 

In Figure 10 (right figures) we show the convergence of our analysis from LO to N 3 LO, which 
demonstrates perturbative stability for the parton distribution functions. The lines rise in the 

6 At 0(a^) only the corrections for Fl(x, Q 2 ) in the region Q 2 ^> nig are known, cf. [44,45]. 
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large x region and drop slightly for low values of x. The shifts from NNLO to N 3 LO are far 
within the present error bands. 

Another way to compare the NNLO fit results consists in forming moments of the distributions 
u v (x, Q 2 ), d v (x, Q 2 ), and u v (x, Q 2 ) — d v (x, Q 2 ). These moments may be directly compared with 
the results of lattice simulations at low pion masses using dynamical quarks in the near future. In 
Table 4 we present the lowest non-trivial moments of these distributions at Q 2 = Q\ in N 3 LO and 
NNLO and compare to the respective moments obtained for the parameterizations [12,13,49]. 
The recent parameterization [13] comes somewhat closer to our present parameterization then 
an earlier one [12] due to an improvement in the x(d(x, Q 2 ) — u(x, Q 2 )) distribution and taking 
into account off-shell effects in the target mass corrections [51]. The errors of the moments 
for [13] are slightly smaller than those in the present analysis. There are still some deviations, 
in particular for d v (x, Q 2 ), but most of the values get very close. 

In Figures 11 and 12 we show the evolution of the valence quark densities xu v (x, Q 2 ) and 
xd v (x,Q 2 ) from Q 2 = 1 GeV 2 to Q 2 = 10 4 GeV 2 in the region x e [1(T 3 , 1] at NNLO with 
correlated errors. To compute the respective error bands error propagation has to be performed 
through the evolution equations. We also compare to other NNLO analyses [12,49]. With rising 
values of Q 2 the distributions flatten at large values of x and rise at low values. Numerically fast 
parameterizations of the distributions xu v (x,Q 2 ) and xd v (x,Q 2 ) and their errors in the region 
Q 2 = [1, 10 6 ] GeV 2 , x e [1(T 9 , 1] are performed and are available from the authors on request. 



6 Aqcd and a s (M|) 



The QCD analysis is performed in NLO, NNLO and N 3 LO. The QCD scale Aq'^ is determined 
together with the parameters of the parton distributions with correlated errors. As an example 
the correlation matrix elements of the NNLO analysis are shown in Table 3. 

In spite of the unknown 4-loop anomalous dimensions one may wonder whether any statement 
can be made on the non-singlet parton distributions and Aqcd at the 4-loop level. The 3-loop 
Wilson coefficients are known [8] and the question arises, which effect has the 4-loop anomalous 
dimension if compared to the Wilson coefficient. The gross behaviour of the 4-loop anomalous 
dimension may be Pade-approximated and we will allow a ± 100% error for it in the analysis. 
After the present analysis was finished, the 2nd moment of the 4-loop anomalous dimension was 
calculated in [10] 
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and can be compared to the Pade-approximant 
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(2)' 



(23) 



For the second moment we found a deviation of about 20% only in the 0(af) term, which is far 
inside our uncertainty margins in the non-singlet case. The variation in Aqcd varying the 4-loop 
anomalous dimension between zero and twice the value of the Pade-approximant amounts to 
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±2 MeV. Thus Aqcd at 0(aj) accuracy is widely determined by the 3-loop Wilson coefficient 
and the behaviour of a s (Q 2 ) to 4-loop order. 

Table 5 summarizes our fit results comparing AJcd 4 and a s (M 2 z ) for the NLO, NNLO and 
N 3 LO analysis. The NNLO value of AqqJ 4 is found to be smaller than the NLO value, while the 
N 3 LO value comes out somewhat higher than the NNLO value at about the same experimental 
error. The pattern of central values from NLO to N 3 LO shows the convergence of the analysis, 
which is illustrated by the decreasing difference of these values. This measure of convergence is 
more significant than varying the scales of the problem, which is less well defined since the range 
of variations is arbitrary. 

We compare the results of the present analysis to results obtained in the literature at NLO 
and NNLO in Table 6. Most of the NLO values for a s (M 2 z ) shown were determined in combined 
singlet- and non-singlet analyses, partly including jet-data in pp scattering. Not in all analyses 
the correlation matrices for the parameters of the fit were presented. The NLO values for 
a s (M z ) are larger than those at NNLO as several analyses, which performed both fits, show. 
The difference of both values, however, is not always the same, most likely due to the type of 
analysis being performed (singlet and non-singlet, non-singlet only, etc.), in which also partly 
different data sets are analyzed. 

As outlined above we perform kinematic cuts in the present analysis to eliminate as widely 
as possible higher twist effects. In this way a significant part of the SLAC data is cut away 
due to the requirement of a minimal value W 2 = 12.5 GeV 2 . Other analyses, cf. [12,18], fit 
phenomenological models for the higher twist contributions. In a precision determination of 
a s (M z ) we do not follow this way since a thorough description of these contributions in terms of 
the light-cone expansion, including the description of the 4-, 6- etc. parton correlation functions, 
the corresponding Wilson coefficients and anomalous dimensions, is not known yet, but would 
be required. 

A significant part of the data in the non-singlet case is due to the BCDMS measurements [2] . 
These data have repeatedly being criticized to deliver a too small value for Aq CD , respectively 
a s (M z ). Let us therefore revisit associated analyses. The NLO BCDMS-fit was confirmed 
in [52] both in the singlet and a separate non-singlet analysis with practically the same value 



for AS' MS , 



A J C d' MS = 198 ± 20 MeV ; a s (M 2 z ) = 0.1096 . (24) 



In a re-analysis of the BCDMS data [53] using somewhat less severe y cuts than applied in the 
present paper the corresponding value 

Ajco 4 '^ = 257 ± 40 MeV, (25) 

was found, which is numerically in accordance with our result. However, in [53] the flavor 
threshold matching in a s (/i 2 ) was performed at the rather high scales 4mg, m 2 = 9 GeV 2 , m 2 = 
80 GeV 2 . Usually the convention is to match at m c and m^. Therefore, the results of [53] cannot 
be compared to other analyses directly. We preformed the NLO analysis with these matching 
conditions and obtained 

Aj c = D 4 ^ = 254±26MeV; a s (M 2 ) = 0.1127 ( 26 ) 
and do not confirm the rather high values of a s (M z ) reported in [53]. 
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Non-singlet QCD analyses were also performed for neutrino data. In [15] the CCFR iron 
data on xF 3 (x,Q 2 ) [54] were analyzed in NLO and NNLO using fixed moments. Likewise a 
NNLO analysis was performed in [14]. Iron data suffer form the EMC effect, Fermi motion and 
other nuclear effects, as re-scattering, for which the QCD scaling violations are not known in 

detail. In [15] rather large values for Aqq,}' ms 

AQ f c = D 4 So = 371±72MeV (27) 
A$gg? LO = 316±51MeV (28) 

are obtained, which are higher than the values obtained in the analyses based on F$' d (x, Q 2 ) 
data, still showing the pattern that the NNLO value is lower than the NLO value. In a quite 
similar analysis in Ref. [14] on the other hand, values close to the results of the present analysis 
are obtained, 

AJcdS = 281±57MeV (29) 
ASt = 255±55MeV. (30) 

In Figure 13 we compare the results of different NLO and NNLO QCD unpolarized and 
polarized DIS analyses for a s (M|) with the results obtained in the present paper and the world 

average for a s (M|) 

a s (Mf) = 0.1189 ±0.0010 , (31) 

which has been determined in [55] recently. In most of the cases the central values lay below 
the world average. Yet most of the NLO results are la compatible with the world average and 
an increasing number of NNLO results departs to lower values. In the case of polarized deep 
inelastic scattering at present only NLO analyses are performed with larger errors than in the 
unpolarized case, which is due the current experimental errors of the polarization asymmetries. 

Nf — 2 MS 

We finally would like to mention results on Aq'^' from lattice calculations with two dy- 
namical quark flavors [61,62] 

Aq^ 2 ' 5 ^ = 245 ± 16(stat.) ± 16(syst.) MeV (32) 
AQ f c = D 2,M5 = 261 ± 17(stat.) ± 26(syst.) MeV . (33) 

These values are very close to the Nf = lattice results with rather small errors. It will be 
interesting to see which values will be obtained for four dynamical quark flavors in the future. 



7 Higher Twist 

The cuts applied in Q 2 and W 2 , Q 2 > 4 GeV 2 and W 2 > 12.5 GeV 2 in the standard analysis 
do widely eliminate higher twist effects from the data sample. Still the region for deep-inelastic 
scattering ranges to lower values of W 2 > 4 GeV 2 keeping the cut in Q 2 . The data in the 
kinematic region 

Q 2 > 4 GeV 2 , 4 < W 2 < 12.5 GeV 2 (34) 
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may thus be used to yield information on power corrections in F%{x, Q 2 ) and F$ (x, Q 2 ). For this 
purpose we extrapolate the QCD fit results obtained for W 2 > 12.5 GeV 2 to the region (34) and 
form the difference between data and theory, applying target mass corrections in addition. The 
empirical higher twist coefficient 

Fr(x, Q 2 ) = O tmc [F? 2 (x, Q 2 )} ■ (l + °q^ 2 ^ ) (35) 

is extracted. Here the operation Otmc[--] denotes taking the target mass corrections of the twist- 
2 contributions to the respective structure function. The coefficients Cut(x, Q 2 ) are determined 
in bins of x and Q 2 and are then averaged over Q 2 . The extracted distributions for Cht(^) 
are depicted in Figures 14 and 15 for the non-singlet case considering scattering off the proton 
and the deuteron target, respectively. The coefficient Cwr(x) grows towards large x. We find 
a similar pattern if compared to the early NLO analysis [63]. Comparing the extractions for 
the case of calculating F% w2 (x, Q 2 ) in NLO, NNLO and N 3 LO a gradual reduction of C U t{x) 
is obtained, i.e. higher twist extractions describing the twist-2 contributions by low order in 
perturbation theory only lead to incorrect results. As seen from Figures 14, 15 Cht(^) is widely 
independent of the target comparing the results for deeply inelastic scattering off protons and 
deuterons, after nuclear wave function corrections were performed. Beyond x ~ 0.9 the errors 
for Cn T (x) become very large. 

8 Conclusions 

We performed a twist-2 QCD analysis of the non-singlet world data up to NNLO and N 3 LO 
and determined the valence quark densities xu v (x,Q 2 ) and xd v (x,Q 2 ) with correlated errors 
in the MS-scheme using conservative kinematic cuts for the data. Parameterizations of these 
parton distribution functions and their errors were derived in a wide range of x and Q 2 as fit 
results at LO, NLO, NNLO, and N 3 LO. In the analysis the QCD scale and the strong 

coupling constant a s (M|), were determined up to N 3 LO. The effect of the 4-loop anomalous 
dimension, if compared to the 3-loop Wilson coefficient and the 4-loop effects in a s (fi 2 ) turn 
out to be rather weak. The corresponding estimates, based on a Pade-approximation compare 
very well with the recently calculated second moment of the 4-loop anomalous dimension. The 
assumption of a ±100% error on the Pade-approximant results into a marginal shift of ±2MeV 
in Aqq^ only. The convergence of the central values for Aqq^ and a s (M|) from NLO to N 3 LO 
shows the stability of the analysis and provides some measure of remaining higher order effects. 
Moments for the valence quark distributions with correlated errors at NNLO and N 3 LO were 
calculated for comparisons with upcoming lattice simulations. Higher twist effects both in F% 
and F 2 d are extracted extrapolating the twist-2 results to the region W 2 < 12.5 GeV 2 . These 
contributions are found to be widely target independent and decrease with rising order in the 
coupling constant for the twist-2 contributions. 
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1 
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Table 1: Number of data points for the non-singlet QCD analysis with their x and Q 2 ranges. In 
the first column are given (in parentheses) the beam momentum in GeV of the the respective data 
set (number), a flag whether the data come from a combined analysis of all beam momenta (comb) 
or whether the data are taken at high momentum transfer (hQ2). The fourth column (F 2 ) contains 
the number of data points according to the cuts: Q 2 > 4 GeV 2 , IT" 2 > 12.5 GeV 2 , x > 0.3 for 
F$ and F 2 rf and x < 0.3 for F^ s . The reduction of the number of data points by the additional 
cuts on the BCDMS data (y > 0.3) and on the NMC data (Q 2 > 8 GeV 2 ) are given in the 5th 
column (F 2 cuts). The 6th column (F 2 HT) contains the number of data points in the range 
4 GeV 2 < W 2 < 12.5 GeV 2 used to fit the higher twist coefficients Cht(x) for the proton and 
deuteron data. In the last column the normalization shifts (see text) are listed. 
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NLO 


NNLO 


N 3 LO 


u v 


a 
b 

P 

7 


0.274 ± 0.027 
3.909 ± 0.040 

6.003 

35.089 


0.291 ± 0.008 
4.013 ± 0.037 

6.227 

35.629 


0.298 ± 0.008 
4.032 ± 0.037 

6.042 

35.492 


d v 


a 
b 

P 

7 


0.461 ± 0.030 
5.683 ± 0.228 

-3.699 

16.491 


0.488 ± 0.033 
5.878 ± 0.239 

-3.639 

16.445 


0.500 ± 0.034 
5.921 ± 0.243 

-3.618 

16.414 


Aqcd 5 MeV 


265 ± 27 


226 ± 25 


234 ± 26 


X 2 /ndf 


484/546 = 0.89 


472/546 = 0.86 


461/546 = 0.84 



Table 2: Parameter values of the NLO, NNLO and N 3 LO non-singlet QCD fit at Q 2 = 4 GeV 2 . The 
values without error have been fixed after a first minimization since the data do not constrain these 
parameters well enough (see text). 



NNLO 


A N f =4 
A QCD 




K v 




bd v 


A(4) 
^QCD 


6.45E-4 










a u v 


9.03E-5 


5.75E-5 








K v 


-3.37E-4 


1.55E-4 


1.40E-3 








1.92E-4 


-8.97E-6 


-4.69E-4 


1.07E-3 




bd v 


9.19E-4 


5.82E-5 


-3.30E-3 


7.21E-3 


5.72E-2 



Table 3: The covariance matrix of the NNLO non-singlet QCD fit at Ql = 4 GeV 2 . 
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/ 


n 


N 3 LO 


NNLO 


MRST04 


A02 


A06 


u v 


2 


0.3006 ± 0.0031 


0.2986 ± 0.0029 


0.285 


0.304 


0.2947 




3 


0.0877 ± 0.0012 


0.0871 ± 0.0011 


0.082 


0.087 


0.0843 




4 


0.0335 ± 0.0006 


0.0333 ± 0.0005 


0.032 


0.033 


0.0319 


d v 


2 


0.1252 ± 0.0027 


0.1239 ± 0.0026 


0.115 


0.120 


0.1129 




3 


0.0318 ± 0.0009 


0.0315 ± 0.0008 


0.028 


0.028 


0.0275 




4 


0.0106 ± 0.0004 


0.0105 ± 0.0004 


0.009 


0.010 


0.0092 




2 


0.1754 ± 0.0041 


0.1747 ± 0.0039 


0.171 


0.184 


0.182 




3 


0.0559 ± 0.0015 


0.0556 ± 0.0014 


0.055 


0.059 


0.057 




4 


0.0229 ± 0.0007 


0.0228 ± 0.0007 


0.022 


0.024 


0.023 



Table 4: Comparison of low order moments at Ql = 4 GeV 2 from our non-singlet N 3 LO and NNLO 
QCD analyses with the NNLO analyses MRST04 [49], A02 [12] and A06 [13] derived from global 
analyses, resp. combined singlet/non-singlet fits. 
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AqFd, MeV 


a s (Ml) 


1NLU 


ZOO ±21 


0.1148 


+0.0019 
-0.0019 


(expt) 


NNLO 


226 ± 25 


0.1134 


+0.0019 
-0.0021 


(expt) 


N 3 LO 


234 ± 26 


0.1141 


+0.0020 
-0.0022 


(expt) 



Table 5: Aj c = 4 and a s (M z ) at NLO, NNLO and N 3 LO. 





a s (M 2 z ) 


expt 


theory 


model 


Ref. 


NT,D 








0.1165 


±0.0065 






M Ql 


mr QTn^ 
ivirto 1 uo 


0.1165 


±0.0020 


±0.0030 




U 11 


A02 


0.1171 


±0.0015 


±0.0033 




[12] 


ZEUS 


0.1166 


±0.0049 




±0.0018 


[56] 


HI 


0.1150 


±0.0017 


±0.0050 


+0.0009 
-0.0005 


[5] 


BCDMS 

1 J V. 1 y 1V± kJ 


0.110 


±0.006 






[21 


GRS 


0.112 








[18] 


BBG 


0.1148 


± 0.0019 








BB (pol) 


0.113 


±0.004 


+0.009 
-0.006 




[58] 


NNLO 






MRST03 


0.1153 


±0.0020 


±0.0030 




[11] 


A02 


0.1143 


±0.0014 


±0.0009 




[12] 


SY01(ep) 


0.1166 


±0.0013 






[14] 


SY01(i/N) 


0.1153 


±0.0063 






[14] 


GRS 


0.111 








[18] 


A06 


0.1128 


± 0.0015 






[13] 


BBG 


0.1134 


+0.0019 
-0.0021 








N 3 LO 






BBG 


0.1141 


+0.0020 
-0.0022 









Table 6: Comparison of a 8 (M z ) values from NLO, NNLO, and N 3 LO QCD analyses. 
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10 Figures 




Figure 1: The MRST parameterization of (d — u) [40] describing the E866 data [39]. Here also 
compared with a parameterization given in [12]. [Courtesy by J. Stirling.] 
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Figure 2: Relative size of the NLO contribution due to cc-production for m c = 1.5 GeV for 
F%{x, Q 2 ) in the valence quark region. 
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Figure 3: Relative size of the NLO contribution due to cc-production for m c = 1.5 GeV for 
F 2 d (x, Q 2 ) in the valence quark region. 
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Figure 4: Relative size of the NLO contribution due to cc-production for m c = 1.5 GeV for 
F 2 NS (x, Q 2 ) = 2[F 2 p (x, Q 2 ) - F$(x, Q 2 )] in the region 0.01 < x < 0.3. 
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Figure 5: The structure function F% as function of Q 2 in intervals of x. Shown are the pure 
QCD fit in NNLO (solid line) and the contributions from target mass corrections TMC (dashed 
line) and higher twist HT (dashed-dotted line). The arrows indicate the regions with W 2 > 
12.5 GeV 2 . The shaded represent the fully correlated la statistical error bands. 
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Figure 6: The structure function F$ as function of Q 2 in intervals of x. Shown are the pure 
QCD fit in NNLO (solid line) and the contributions from target mass corrections TMC (dashed 
line) and higher twist HT (dashed-dotted line). The arrows indicate the regions with W 2 > 
12.5 GeV 2 . The shaded represent the fully correlated la statistical error bands. 
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Figure 7: The structure function F 2 as function of Q 2 in intervals of x. Shown is the pure QCD 
fit in NNLO (solid lines) with its fully correlated la statistical error band (shaded areas). 
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Figure 8: The structure function F 2 NS as function of Q 2 in intervals of x without data points. 
Shown are the pure QCD fit in NNLO (solid lines) with its fully correlated la statistical error 
band (shaded areas) and the contribution from target mass corrections TMC (dashed line). 
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Figure 9: The parton densities xu v and xd v at the input scale Q% = 4.0 GeV 2 (solid line) 
compared to results obtained from NNLO analyses by MRST (dashed-dotted line) [49] and A06 
(dashed line) [13] (left panels) and from NLO analyses by CTEQ6 (dashed-dotted line) [19], 
GRV (dashed line) [50] and HI (long dashed line) [5] (right panels). The shaded areas represent 
the fully correlated la statistical error bands. 
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Figure 10: Left panels: The parton densities xu v and xd v at the input scale Ql = 4.0 GeV 2 (solid 
line) compared to results obtained from NNLO analyses by MRST (dashed-dotted line) [49] and 
A06 (dashed line) [13]. The shaded areas represent the fully correlated la statistical error bands. 
Right panels: Comparison of the same parton densities at different orders in QCD as resulting 
from the present analysis. 
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Figure 11: The parton density xu v at NNLO evolved up to Q 2 = 10,000 GeV 2 (solid lines) 
compared to results obtained by MRST (dashed-dotted line) [49] and A02 (dashed line) [12]. 
The shaded areas represent the fully correlated la statistical error bands. 
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Figure 12: The parton density xd v at NNLO evolved up to Q 2 = 10, 000 GeV 2 (solid lines) 
compared to results obtained by MRST (dashed-dotted line) [49] and A02 (dashed line) [12]. 
The shaded areas represent the fully correlated la statistical error bands. 
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Figure 13: Comparison of different NLO, NNLO and N 3 LO measurements of the strong coupling 
constant a s (M§) in unpolarized and polarized deeply inealstic scattering. CTEQ6 [19], MRST03 
[11], A02 [12], ZEUS [56], HI [5], BCDMS [2], BBG: present analysis, SY01 [14], A06 [13], BB [58], 
SMC [59], ABFR [60]. 
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Figure 14: The higher twist coefficient Cht{x) for the proton data as function of x treating 
twist-2 contributions in NLO, NNLO and N 3 LO. 
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Figure 15: The higher twist coefficient Cht{x) for the deuteron data as function of x treating 
the twist-2 contributions in NLO, NNLO and N 3 LO. 
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